Spontaneous and induced dynamic correlations in glass-formers II: 
Model calculations and comparison to numerical simulations 
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We study in detail the predictions of various theoretical approaches, in particular mode-coupling 
theory (MCT) and kinetically constrained models (KCMs), concerning the time, temperature, and 
wavevector dependence of multi-point correlation functions that quantify the strength of both in- 
duced and spontaneous dynamical fluctuations. We also discuss the precise predictions of MCT 
concerning the statistical ensemble and microscopic dynamics dependence of these multi-point cor- 
relation functions. These predictions are compared to simulations of model fragile and strong 
glass-forming liquids. Overall, MCT fares quite well in the fragile case, in particular explaining the 
observed crucial role of the statistical ensemble and microscopic dynamics, while MCT predictions 
do not seem to hold in the strong case. KCMs provide a simplified framework for understanding how 
these multi-point correlation functions may encode dynamic correlations in glassy materials. How- 
ever, our analysis highlights important unresolved questions concerning the application of KCMs to 
supercooled liquids. 
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I. INTRODUCTION 

Dynamic heterogeneity is a well established feature of 
the behavior of a diverse class of systems close to their 
glass transition temperatures [l|, S H, 0, H, @|. Given 
the relatively recent realization of the importance of dy- 
namic heterogeneity, it is not surprising that the system- 
atic characterization of such spatio-temporal behavior, 
and the Icngthscalcs associated with it, is far from com- 
plete. Much recent effort has been expended to devise 
correlation functions that effectively and quantitati vely 
probe dynamic heterogeneity 0, i, H [13, [II!, El, [H, [li- 
The theoretical understanding of the behavior of such 
correlation functions is still in its infancy. This is to 
be contrasted with our relatively mature understanding 
of bulk structure and dynamics in supercooled liquids 
as measured by simple, low order correlation functions 
(such as intermediate scattering functions) that can only 
indirectly hint at dynamic heterogeneity P, ■ 

In the first paper of this series, denoted in the follow- 
ing as I [l6t , we set out to provide a general understand- 
ing of the behavior of a particular class of multi-point 
correlation functions that encode information concerning 
the growing dynamical lengthscale in supercooled liq- 
uids. To set the stage for the present work we briefly 
recall some definitions and results obtained in I. Let 
/(r,t) = o(r, <:)o(r, 0) be the instantaneous value of a 
local two-time correlator at position r and time t, and 
[/(i)]r = J (i''r/(r,t) its spatial average over a large 
but finite volume V. The thermal average {[f{t)]r) is 
a standard two-time correlator, such as the intermedi- 
ate scattering function when the observable o(r, t) is the 



excess density p{v,t) — pQ. A previously defined multi- 
point susceptibility is the following four-point dynamic 
susceptibility, 

Xi{t)^N{[d.m]l)=p f d''r{5f{r,t)5mt)), (1) 



where we introduced the notation 6X = X — {X) for the 
fluctuations of the observable X. From Eq. ([T|), we see 
that Xii^) quantifies the strength of spontaneous fluc- 
tuations of the dynamical behavior in supercooled liq- 
uids by their variance. As shown by the last term in 
Eq. ([T|) fluctuations become larger if this dynamic hetero- 
geneity becomes increasingly spatially correlated. Since 
Xi{t) is the volume integral of the four-point correlator 
S4{r,t) ~ ((5/(r, t)(5/(0, t)) (or, alternatively, in Fourier 
space, X4{t) = liniq^o 'S'4(q, f)), it is directly related 
to the number of correlated particles, Xii^) ^ i^,/'^)'^^ j 
where ^ is the dynamic correlation length, a a molecu- 
lar lengthscale, and d/ is related to the possibly fractal 
geometry of the dynamic heterogeneity. The direct link 
between Xiit) and the lengthscale of dynamic heterogene- 
ity ^ explains the intensity of the present experimental 
effort dedicated to its measurement [TTb J^, 2^ ■ 

Recognizing that spontaneous fluctuations are in gen- 
eral hard to access experimentally, we have suggested |16l . 
[Tgl to measure instead the response of the averaged two- 
time dynamical correlators to an infinitesimal perturbing 
field. 



diifjt)] 
dx 



(2) 



In particular we have dedicated much effort to the cases 
where x is either the temperature, a; = T, or the density. 
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X = p, focusing therefore on Xt(0 ^^'^ Xpi^)- In R-sf. [T^ 
it was argued that Eq. ^ defines an experimentally ac- 
cessible multi-point dynamic susceptibility which is a rel- 
evant alternative to Xii^)- There are two important ar- 
guments to support this claim, which we only summarize 
for XT{t), but they directly carry out to Xp(i). The first 
one is that, for a classical fluid evolving via Newton's 
equations at constant number of particles, N, volume V, 
and energy, E, the following fluctuation-dissipation the- 
orem holds: 

kBr'xT{t)^V{[Sf{t)]r[dh{t)]r) - J d''r{5f{r,t)6h{0,0)), 

(3) 

where [e(t)]r = V ^ J <i'^re(r, i) is the instantaneous 
value of the energy density, and fcs the Boltzmann con- 
stant. The similarity between Eqs. H]) and ^ is strik- 
ing. The new susceptibility xrit) quantifies the strength 
of correlations between dynamic fluctuations and energy 
fluctuations. As shown by the last term in Eq. ([3]) XT{t) 
becomes larger if dynamical and energy fluctuations be- 
come increasingly spatially correlated. Since XT{t) is 
proportional to the volume integral of the three-point 
correlator ^^(r, i) — {Sf{r,t)Sh{0,0)) (or, alternatively, 
XT{t) = hmq^o S'T(q, t)), it is also directly related to a 
correlation volume, which makes it an equally appealing 
quantity. A second argument establishing the relevance 
of XT{t) is the fact that xt(0 and X4{t) can be related 
by the following inequality: 

X4{t) > —T\Ut), (4) 

CP 

where cp — V{[Sh{t)]l)/{pT'^) is the constant volume 
specific heat expressed in units of ks- The result (j4|) 
can be understood by formal consideration about statis- 
tical ensembles (see I), or more simply by noting that the 
relation ^ stems from the fact that the (squared) cross- 
correlation between two observables (encoded in xt(0) 
cannot be larger than the product of their variances (en- 
coded in Xi(t) and Cp). 

In I, we focused on the thermodynamic ensemble de- 
pendence and the dependence on the microscopic dy- 
namics. Using general theoretical arguments, we gave 
qualitative and quantitative guidelines for these depen- 
dences. The ensemble variability of global multi-point 
indicators of dynamical heterogeneity (corresponding to 
fluctuations of intensive dynamical correlators) is not sur- 
prising, given what is already understood about the en- 
semble dependence of sirtipler susceptibilities near stan- 
dard critical points [2ll . [2^ . Importantly, this ensem- 
ble dependence allows for the derivation of the rigorous 
bound ^ on Xii't) that is potentially useful for providing 
a simple experimental estimate of the lengthscale asso- 
ciated with dynamical heterogeneity near Tg. That this 
bound becomes a good approximation for Xii^) above 
Tg was checked in simulations of both strong and fragile 
glass forming liquids in I. The predicted dependence on 
the underlying nature of the dynamics is perhaps more 



surprising, especially in light of the fact that simulations 
of simple dynamical correlation functions show no such 
non-trivial dependence [I^IU, H^. Again, in I we have 
confirmed this striking prediction by atomistic simula- 
tions. 

Having outlined some generic properties of a class of 
multi-point indicators of dynamical heterogeneity, and 
confirmed these basic predictions in I, we now turn to 
the information contained in specific theories of glassy 
dynamics. In particular, we address in this paper various 
properties of these susceptibilities from the standpoint 
of simple mean- field spin-glass models [2^, the mode- 
couphng theory (MCT) of supercooled liquids [27| . and 
kinetically constrained models (KCMs) [23 • Our choice 
of theoretical models is natural: to our knowledge, only 
MCT and KCMs offer a detailed theoretical description 
of dynamic heterogeneity in supercooled liquids. We aim 
to confront these theories with the general theoretical 
properties outlined in I, as well as with simulations of 
atomistic glass-forming systems. The outcome of this 
exercise will be a greater understanding of the successes 
and failures of these theories, which will lead us to for- 
mulate a number of questions related to the comparison 
of these models with the expectations outlined in I. 

This paper is organized as follows. In Sec. [Til we 
present the results predicted by MCT. This includes gen- 
eral scaling behavior, as well as the dynamics and ensem- 
ble dependence as derived via a field-theoretical approach 
to MCT. Within this approach we can show in particu- 
lar that strong ensemble and dynamics dependence of dy- 
namic fluctuations arises, while no such dependence is ex- 
pected for averaged quantities [2^ . This section discusses 
the wavelength dependence of Xii^) for liquid state MCT, 
while the relationship between Xii^) and xrit) within 
p-spin models for which MCT is exact is performed in 
Appendix A. In Sec. Illll we turn to the ensemble and dy- 
namics dependence of Xiit) and Xrit) in KCMs. Here, 
we discuss different models with varying degrees of co- 
operativity. Interesting unresolved questions, concerning 
the relevance of KCMs to model molecular glasses, are 
outlined in this section. In Sec. IIVI the prediction of 
these various models are compared to atomistic simula- 
tions. In Sec.[Vl we conclude and we detail the successes 
and failures of the theoretical models in light of the com- 
parison with simulations, and give a summary of these 
comparisons in Table I. 



II. MODE-COUPLING THEORY OF 
DYNAMICAL FLUCTUATIONS 

A. MCT and dynamic fluctuations 

Because it starts from a microscopic description of su- 
percooled liquids and ends up with a complete descrip- 
tion of its dynamics, MCT is a powerful tool for the in- 
terpretation and prediction of the qualitative and quanti- 
tative behavior of slow dynamics in glass-forming liquids 
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and colloids, at least not too close to the glass transi- 
tion [l^l- The MCT transition is usually described as a 
small scale phenomenon, the self-consistent blocking of 



the particles in their local cages [27[. This is surpris- 
ing since on general grounds a diverging relaxation time 
is expected to arise from processes involving an infinite 
number of particles (leaving aside the case of quenched 
obstacles) (30[. Actually, the cage mechanism requires 
some kind of correlation in space: in order to be blocked 
by one's neighbors, the neighbors themselves must be 
blocked by their neighbors and so on until a certain scale 
that, intuitively, sets the relaxation timescale of the sys- 
tem. The fact that within MCT "cages" are correlated 
objects [m will in fact become clear below. 

This "local-cage" point of view was challenged in the 
context of mean-field disordered systems by Franz and 
Parisi [ll|; see also [H, [11] for early results. For these 
models the dynamical equations for correlators are for- 
mally equivalent to the schematic version of the MCT 
equations. Franz and Parisi [ll| argued that a dynami- 
cal susceptibility similar to X4(0 i^i these models has a 
diverging peak at the dynamical mode-coupling transi- 
tion. The Franz-Parisi susceptibility is further discussed 
in Appendix A. Although a lengthscale cannot be de- 
fined in mean- field models, a diverging susceptibility is 
the usual mean-field symptom for a diverging lengthscale 
in finite dimensions. More recently, two of the authors 
(BB) [U, using a field-theoretical approach to MCT, 
clearly showed the existence of a diverging length within 
MCT and analyzed the critical properties of dynamical 
fluctuations. In that work the role of conserved quanti- 
ties, emphasized in I, was overlooked. As we show in the 
following, BB's results for X4{t) are correct either for dy- 
namics without any conserved variables (as is the case for 
disordered p-spin systems with Langevin dynamics), or 
in ensembles where all conserved variables are fixed, i.e. 
NVE for Newtonian dynamics and NVT for Brownian 
or Monte-Carlo dynamics. 

When there are conserved variables the four-point 
correlation function 6*4 (q,t) can be decomposed in two 
terms, in agreement with the general considerations of I. 
These two terms reflect different physical contributions 
for q = 0: one is the contribution in the ensemble where 
all conserved variables are strictly fixed, and the sec- 
ond arises from the fluctuations of dynamically conserved 
variables that feed back into the dynamical correlations. 
The second term (for q = 0) is therefore absent in an en- 
semble where these variables are fixed. This latter term is 
the one that yields a lower bound for limq^o 'S'4(q, i), as 
expressed in Eq. The bound involves the derivative 
Xx{t) defined in Eq. ([S]), where x is a conserved variable. 
For example x = p for hard-spheres, where the density is 
a conserved quantity both for Brownian and Newtonian 
dynamics, or x = H (or £'), the enthalpy (or the energy), 
in cases where temperature is the relevant control param- 
eter. One can of course also focus on dynamical responses 
with respect to thermodynamic control parameters such 
as the pressure or the temperature. One formulation is 



related to the other via a trivial thermodynamic change 
of variables and the chain rule. In the following, for sim- 
plicity, we will always focus on the derivative with respect 
to conserved degrees of freedom. 

In the next subsections we shall uncover the critical 
properties of the dynamical fluctuations and dynamical 
responses discussed above, and obtain and analyze quan- 
titative predictions for dynamical responses within MCT. 
We numerically confirm these results within the p-spin 
model in [Appendix A[ 



B. Dynamic scaling and critical behavior 

In the following, using the field-theoretical framework 
developed in I, we obtain the critical behavior of dynam- 
ical fluctuations close to the MCT transition. We focus 
in particular on X4(t), S'4(q, i), and Xx{t)- 



1. Ladder diagrams within MCT 

Different derivations of MCT follow a common strat- 
egy: write down exact or phenomenological stochastic 
equations for the evolution of the slow conserved degrees 
of freedom and then use a self-consistent one-loop ap- 
proximation to close the equations. For instance, in the 
case of Brownian dynamics the only conserved quantity 
is the density, and the so-called Dean-Kawasaki equa- 
tion [H, 111] has been analyzed (see Refs. [13, [38[. l39l] 
for a discussion of the different flcld-theories) . Field- 
theories are obtained through the Martin-Siggia-Rose- 
dcDominicis-Janssen method, where one first introduces 
response fields enforcing the correct time evolution and 
then averages over the thermal noise [il] . 

The direct derivation of MCT equations starting from 
field-theory is difficult and different approaches have been 
pursued [38[. It is still unclear how to obtain in a con- 
sistent way the standard MCT equations derived by the 
Mori-Zwanzig formalism [53, [H, |4l|. Indeed, if time- 
reversal symmetry is preserved, one-loop self-consistent 
equations are not the standard MCT equations, but have 
similar qualitative properties [13]. They lead in particu- 
lar to the same critical behavior of the correlators. This 
issue is not relevant here because we focus on qualitative 
properties of dynamic fluctuations which depend only on 
the critical properties of the MCT transition. 

The starting point for describing dynamic fluctuations 
within field theory is the Legendre functional [H, [Z^, 
[43! r(5'a,Ga,(,) (here and in the following we use the 
notations introduced in I): 

r(*a,Ga,6) = -iTrlogG+iTrG^ilG+**]-$2P/(*,G), 

where '^2Pi{^ ai Gafi) is the sum of all two-particle irre- 
ducible Feynman diagrams (namely those that cannot be 
decomposed in two disjoint pieces by cutting two lines) 
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FIG. 1: Three diagrams approximating $2P/(^a, Ga.b) within 
the MOT approximation. The resulting expression of the self- 
energy is also shown. Lines are full propagators, and dots are 
conserved variables. 
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FIG. 2: Diagrammatic expression of (5^'I>2p//(5'I'i(5G2,3 and 
(5^$2P7/5Gi,25G3,4 within MOT. 



constructed with the vertices of the theory, using the full 
propagator G as the lines and 5* as the sources (Go is 
the bare propagator) [iO, H^. The first derivative of 
r(^'a,Ga.h) leads to the self-consistent equations for the 
order parameters themselves (including G's), whereas the 
second derivatives lead to the equation for the fluctua- 
tions of the order parameters. 

All field-theoretic derivations of MCT consist of one- 
loop self-consistent equations for the dynamical structure 
factor. At the level of the functional this corresponds to 
an approximation of ^2Pi{'^ a,Ga,b) in which only the 
first three diagrams of Fig. [T] are considered. They are 
constructed from a three-leg vertex that is present in all 
field-theories of dense liquids. The black dots represent 
the 5'^ attached as sources and the lines the full prop- 
agators of the theory. The corresponding expression of 
the self-energy E = S^/5G is also shown. Note that the 
second diagram is not present in the usual expression 
for the self-energy used in the MCT equations because 
the solution of the self-consistent equation for leads to 
S'i' = 0. As discussed in I, the reason is that the average 
value of the response fields is zero and the bare values of 
the physical slow fields are not corrected at any order of 
the self-consistent expansion because they correspond to 
conserved variables, whose average value is not fixed by 
the dynamics but through the initial conditions. 

However, when the matrix of second derivatives of T is 
considered, it is important to keep the second self-energy 
diagram because it gives a contribution that cannot be 
neglected. The matrix of second derivatives reads: 



5^r 



l5Gl,2(5G3,4 
S^lSG2.3 " 



/" i — 1 /- ' — 1 



='2P/ 



(SGi,2<5G, 



3,4 



^2PI 

" 5^i5G2,3'' 

(Go')l,2. 



(6) 



As in I, we denote these operators respectively as A, 
B, and G. The diagrammatic expressions for the sec- 



ond derivatives of $2P/ are shown in Fig. [51 Note 
that we show only the contributions that are non-zero 
when evaluated for the average quantities (in particular 
for J* = 0). 

All dynamic fluctuations can then be expressed in 
terms of A^^,B,C (see I). In particular, the four-point 
fluctuations 

( { (X, t)4,t, (X' ,t')-^>a (X, t) (X , i') } 

X {My,s)My',s')-^c{y,s)^diy',s')})^, (7) 

are given by 

A-^ + {A-^B){C - B^A-^B}-^{A-^B)^, (8) 

evaluated at the matrix element [a, 6, x, x',f,f'; 
c,d,y,y' , s, s']. The explicit expression for A~^,B,C 
makes it clear that within MCT the critical properties 
of dynamical fluctuations come only from A~^. Indeed, 
G is just the inverse of the bare propagator, whereas B 
is the bare vertex. These quantities have no critical be- 
havior at the MCT transition. Instead, using the general 
results of I and the MCT expression of 5'^^2Pi /SG6G, 
one finds that A~^ corresponds to the sum of n-ladder 
diagrams shown in Fig. [3] As shown in Ref. [U, the 
resummation of these diagrams indeed leads to a critical 
contribution at the MCT transition. 

In particular, they give a contribution to the four 
point function (Sp-kAt)Spk,+o(0)Sp^k,(t)5pk,,^a(0)) 
that scales as Isil, [sS : 



1 



/e + q 
1 



t 

^ ^ ■ /e Tg 

t 

9a 



t (X Tp, 



t (X Ta. 



(9) 



Note that here and in the following we use the standard 
MCT notation [13, El. In particular e = \xc - x\/xc 
is the reduced distance from the critical point at a; = 
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5^ 
5G 



GG= 



FIG. 3: Expression for within MCT. It consists in 

a sum of n-ladders constructed from the elementary block 
(STi / 5G)GG shown in the second line, see Eq. (6-a). 



Xc- The function i^q^ / ^/'^p) behaves as {t/rp)"' and 
it/rpf' for small and large value of {t/rp), respectively. 
Furthermore the function (t/Ta) behaves as {t/Taf' for 
small values of {t/Ta), see also Eqs. ^Bl [TT]) below. 

Below we show that the critical behavior that emerges 
from ladder diagrams underlies all of the critical proper- 
ties encoded in the dynamical fluctuations and responses, 
as discussed in general terms in I. 



2. Dynamical responses 

Dynamical response functions, defined in Eq. ([2]), are 
particularly interesting because they provide, through in- 
equalities like Eq. (HI) , an estimate of the relevant dynam- 
ical fluctuations, and because they are related to three- 
point dynamical correlations, Eq. ([3]). 

An exact expression for dynamical response functions 
can be derived noting that G is obtained by setting 
dT/dG = 0. Differentiating this relation with respect to 
a conserved variable \1/ one can easily derive the relation 
(see Eq. (60) of I): 



dG 



n -1 



dGdG 



dGd^ 



(10) 



which is represented in Fig. [4] using MCT diagrams. 

Since the derivative is taken with respect to the average 
value of one of the conserved degrees of freedom the 
wavevector entering into the ladder diagrams is zero. As 
a consequence, the scaling of dynamical response func- 
tions is given by Eq. setting q = 0: 



B{q 



0) ft 
— ,9/3 — 

\Tf3 



B(q - 0) 



t 



t (X Tp, 



t (X Tr, 



(11) 



where we have dropped the first argument of gp, equal to 
zero here, and B{q = 0) reminds us that there is an ad- 
ditional contribution from the vertex B. We show below. 



FIG. 4: The dynamical response, obtained from Eq. (|10p by 
noting that the inversion oi A = dco^ involves resumming 
ladders that close on flc^F. 



see Eqs. JTH [151 Ull [H]), that these results can alter- 
natively be obtained analytically using standard MCT 
results. However, the field-theoretical derivation shows 
more clearly the role of the ladder diagrams, and is 
crucial to understand the relationship between dynam- 
ical response and dynamic fluctuations. We note that 
from the diagrammatic expression for X'b & clear rela- 
tionship with three-point dynamical correlators appears. 
The diagrammatic expression of the correlation between 
the fluctuation of the dynamical structure factor and 
the fluctuation of a conserved variables reads (see I): 
~A^^B{G — B^^^^i?}"^, which contains the same dia- 
grams as X* ! with a propagator attached at the end [9l] | . 

In order to probe the spatial dependence of dynami- 
cal fluctuations related to ladder diagrams, zero-wavector 
response functions such as x*(*) ^-^e not sufficient, and 
one should consider instead the response to a spatially 
modulated external field [3l|. It was recently proved 
in Ref. (sij that such a g-dependent dynamical response 
function has the same scaling as the one anticipated from 
ladder diagrams in Eq. 



'ence of 



Ensemble and microscopic dynamics 
fluctuations 



In the following we illustrate, within MCT, the de- 
pendence on statistical ensembles and on microscopic 
dynamics of dynamical fiuctuations which we have dis- 
cussed in full generality in I. In particular, one finds that 
although S'4(q, <) and its g — ^ limit are ensemble in- 
dependent quantities, Xiit) = Sii'i ~ OjO does depend 
on the ensemble and on microscopic dynamics. This fact 
refiects the subtle nature of global fiuctuations when the 
thermodynamic limit is taken p^ . 

Applying the general theory developed in I one finds 
that 5*4 (q, t) is given by the ladder diagrams in Fig. [3] 
plus "squared ladders" as shown in Fig. [5] The ladder 
diagrams shown in Fig. [5] are joined by a propagator at 
wavector q. In ensembles where all conserved degrees of 
freedom are fixed, e.g. Newtonian dynamics in the NVE 
ensemble or Brownian dynamics in the NVT ensemble, 
the propagator evaluated at g = vanishes, because con- 
served quantities do not fluctuate on the scale of the sys- 
tem size (and all propagators related to response fields 
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4- Behavior o/S4(q, f) and upper critical dimension 



n 
m 



FIG. 5: Representation in term of diagrams of {A ^ B){C — 
B''A-^B}-^{A-^B)\ see Eq. within MCT. It corre- 

sponds, roughly speaking, to 'squaring' the ladders of Fig. 



are zero because they are proportional to q at small q). 
Therefore, in these ensembles, simple ladder diagrams 
provide the sole contribution to Xiit) within MCT, caus- 
ing Xi{t) to scale as in Eq. ([9]) evaluated at g = 0. This 
is also true for p-spin models, see Appendix A. 

Instead, in ensembles where at least one conserved de- 
gree of freedom is allowed to fluctuate, e.g. the NVT or 
NPT ensembles for Newtonian dynamics, or the NPT 
ensemble for Brownian dynamics, the propagator joining 
the ladders in Fig. [5] does not vanish, and contributes 
only to non-critical prefactors. In this case the diagrams 
corresponding to "squared ladder" diagrams dominate, 
at least close enough to the transition. Note that their 
overall scale might however be small, for example if the 
compressibility or the specific heat are large or if the dis- 
tance to the critical point becomes large. They lead to a 
modified critical behavior for X4(^) within MCT, reading: 



.9/3 , 



■9a 



t OC ra, 

t oc Tq,, 



(12) 



where cja oc and gp oc g^. 

These results provide non-trivial relationship between 
dynamical fluctuations in different ensembles and differ- 
ent microscopic dynamics. Although already suggested 
by the general theory developed in I they become sharp 
statements within MCT. For instance, the predicted dy- 
namic scaling of Xi{t) in the NVT ensemble for New- 
tonian dynamics is that of the "squared ladder" dia- 
grams of Eq. whereas the predicted scaling of Xi{i) 
in the NVT ensemble for Brownian dynamics is that 
of Eq. (fTT|) . which coincides with the one expected for 
Newtonian dynamics in the NVE ensemble. Note that 
the critical mechanism underlying the dynamical fluctu- 
ations is the same for all ensembles and dynamics and is 
uniquely encoded in ladder diagrams. Indeed there is a 
unique lengthscale, diverging in the same way for all mi- 
croscopic dynamics, which underlies the critical behavior. 
However, the coupling to conserved degrees of freedom 
may produce a large amplification of global fluctuations. 



The behavior of S'4(q, t) for q ^ Qis apparently simpler 
because it does not depend on the statistical ensemble. 
Furthermore, all types of diagrams are present in its ex- 
pression so that its qualitative critical behavior is inde- 
pendent of the microscopic dynamics, provided that at 
least density is locally conserved. However, since S'4(q, t) 
contains the two terms discussed previously (ladders and 
squared ladders) a crossover behavior might be expected. 
Although the squared ladders should dominate very close 
to Tc they might become sub-dominant far from the crit- 
ical point. Therefore one should be very cautious when 
comparing the present MCT predictions to the behav- 
ior of real liquids where the mode-coupling singularity Tc 
is replaced by a smooth crossover towards an activated 
regime. In order to judge the relative importance of the 
two terms (ladder and squared ladder) , one may focus on 
their q = {) value, which corresponds to Xi^^ for ladders 
and to kBT'^Xrl cv for squared ladders. For example, in 
the case of the LJ mixture studied in I, the latter term 
becomes dominant only close to the transition T ~ 0.47. 
As a consequence, for higher temperatures, the contribu- 
tion of the squared ladders can be neglected and Si{q, t) 
will have the behavior presented in Eq. ([5]). A similar 
crossover is expected for Xi{t)i confirmed numerically 
in I. The important difference with Si^{(\,t) is that it is 
possible, at least in numerical simulations, to disentangle 
the different contributions to Xii^) by working in differ- 
ent ensembles. 

Although four-point correlators were originally hoped 
to be suited to quantify precisely dynamical hetero- 
geneities in glass-formers, our results show that, al- 
though containing useful information on dynamical het- 
erogeneities, they mix it with other, less interesting phys- 
ical effects. This is a further motivation to study dy- 
namical response to spatially modulated fields introduced 
in (sH . for which only the simple ladder diagrams con- 
tribute and which, therefore, allows one to obtain clearer 
and more direct information on dynamic correlations. In 
future work, it would be extremely interesting to com- 
pare this response function computed within MCT [3l| 
to its direct numerical evaluation in a simulated liquid. 

Finally we note that, within MCT, the scaling of dy- 
namic fluctuations in the ensemble where all conserved 
variables fluctuate is different from the one predicted by 
BB [131 • This implies that the upper critical dimension 
of the theory, found to be dc = 6 in [s^] , has to be re- 
vised accordingly. Focusing on the /5-regime, the fluctu- 
ations of the non-crgodicity parameter in a region of size 
^ £-1/4 grow as 6q oc ^4-d/2^ where d is the space di- 
mension. Imposing, in the spirit of a Ginzburg criterion, 
that 5q must be much smaller than the critical behavior 
of the order parameter, i.e. qc — q ^ v^, one finds that 
fluctuations become dominant below the upper critical 
dimension = 8. In [i^l it is shown that this result 
can be obtained from diagrammatic considerations: be- 
low dc = 8, corrections to MCT are found to diverge in 
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the infrared regime. 

C. The fc-dependence of dynamical fluctuations 
within MCT 

Several different definitions of Y4(t) have been em- 
ployed in the literature [H, M, M, M, M, S H [13, [H] . 
Regardless of definition, since two-point density fluctu- 
ations must depend on the dynamically probed length- 
scale k^^), the detailed behavior of '^iU also 
depend on this lengthscale, see e.g. Physically, the 
dependence of the fluctuations on lengthscale reflects the 
coupling or sensitivity of cooperative motion to behavior 
on the scale of the measured two-point fluctuations. For 
example, the expectation that high-frequency phonons 
do not couple strongly to the large lengthscale dynamic 
heterogeneity is reflected by the fact that a Xii^) that 
focuses on short lengthscales associated with vibrations 
cannot exhibit the sizable normalized peak values that 
are connected to large cooperative lengths (but see the 
discussion in [13] )■ Only at a critical point would one 
expect all modes to couple in such a way that the be- 
havior of Xiit) would exhibit truly universal properties. 
Since this issue cannot be discussed within p-spin models 
which contain no lengthscale (see [Appendix A| , we turn 
instead in this section to liquid state MCT which con- 
tains the complete wavevector dependences of dynamic 
bmctions. 

The dependence of Xii^) 

function of lengthscale 
was first discussed by Glotzer and coworkers [46|, who 
used the definition 

X4W = ^ [{Q'it)) {Q{t)r] , (13) 

where N, V and (3 are the number of particles, the 
volume and the inverse temperature, respectively, and 
Qit) = E.y^^(|r.(0) - r,m and wi\r, ~ r^]) = 1 for 
|ri — r2| < a and is zero otherwise. Here a is a cutoff 
parameter. Lacevic et al. explicitly showed that within 
this definition of X4(i), the value of a that maximizes the 
peak height is close to the global Dcbye- Waller ampli- 
tude of the mean-square displacement [11]. For values 
of the cutoff that are larger or smaller, the absolute am- 
plitude of X4{t) decreases. It can be noted in this work 
that the shape of Xii^) is sensitive to the value of a as 
well, although no systematic study of this dependence 
was investigated. 

Using a different definition of XiW^ namely that de- 
fined by Eq. ([50)1 below, Dauchot et al. noted that for 
a weakly sheared granular system, the slope of XiW in- 
creases as the wavevector decreases [l^ . This particular 
dependence has been studied in more detail in a recent 
work [ist . From both molecular dynamics simulation and 
the direct analysis of a class of kinetic facilitated models. 
Chandler et al. have detailed the lengthscale dependence 
of a variety of definitions of X4(t), and have argued that 
a generic feature of this dependence is that the growth of 



X4{t) to its peak becomes significantly more rapid as the 
intrinsic lengthscale increases. It was also argued [1^ 
that this result is inconsistent with the predictions of 
mode-coupling theory outlined in Ref. [l4j . Since the k- 
dependent Xii^) in facilitated models has been discussed 
in detail in Ref. [i^ , we focus below on the predictions of 
mode-coupling theory. In particular, wc show that, de- 
spite statements to the contrary, mode-coupling theory 
is at least in qualitative accord with the behavior found 
from computer simulation, as detailed in Ref. [48j . 

An important aspect of the physical content of the fc- 
dependence of Xii^) is embodied in the consideration of 
the distinction between /3- and a-relaxation, and the im- 
plication that this distinction holds for the lengthscales 
of dynamic heterogeneity. At a given density and tem- 
perature, a well-defined plateau in the two-point den- 
sity correlator for wavevcctors near the first diffraction 
peak in the static structure factor S{k) will saturate as 
k is decreased. Eventually, as k is decreased further the 
hydrodynamic regime is reached, where the local coop- 
erative processes associated with dynamic heterogeneity 
are averaged out. At a fixed distance from the dynami- 
cal transition temperature Tc, as fc is decreased first the 
/3-relaxation window decreases in duration [i^. Con- 
comitantly, the stretching exponent of the a-relaxation 
increases continuously. Eventually, the crossover is com- 
plete when the /3- window is no longer observable, and the 
stretching exponent saturates at unity, signifying long- 
time hydrodynamic behavior. 

The theoretical considerations made in Ref. [U are 
based on the asymptotic predictions of mode-coupling 
theory. Arbitrarily close to Tc and in systems for which 
Tc is not avoided, the predictions of Ref. [3] are nearly 
universal in the sense that the effects mentioned above 
are only seen in the strict fc — > limit. By consider- 
ing the fc-dependence of induced susceptibility Xx(t) for 
a fixed, finite distance from Tc one gains a qualitative 
understanding of how the universal features expected for 
k near the first diffraction peak of S{k) are changed as k 
decreases. 

As discussed in the previous sections dynamical fluc- 
tuations encoded in ladder diagrams are visible either in 
Xiit) or in the dynamical response. However, dynamical 
responses are accessible to direct quantitative numerical 
computations that are an essential tool in order to dis- 
cuss the crossover issues discussed above. Let us now, 
for completeness and clarity, rederive the results for dy- 
namical responses using only standard MCT results [23| 
being particularly careful about the fc-dependence. 

In the a-regime close to the transition, but still in 
the liquid phase, the dynamical structure factor scales 
like F{k,t) ~ fai't/'Tai^)) where Tq. ~ e~'' and 7 — 
l/2a+l/26. Thus, we find that in the a-rcgime Xx{k, t) = 
dF{k, t) I dx reads: 

xAk,t)^-gl{t/Tc.), 5^(x) = --^. (14) 
Wc temporarily change notation to emphasize the fc- 
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dependence, namely, we promote Xxit) to Xxik, t). In the 
/3-regime and close to the transition F{k, t) ~ S{k)q{k) + 



S{k)h(k)^/lfp{t/Ti3) where Tp 



q{k) is the non- 



ergodic parameter and h{k) is the critical amplitude. 
Thus, we find that in the /3-regime Xx{k,t) reads: 



(15) 



Analyzing the /3-regime with a large but not diverging 
time and matching the a-regime with the /3-regime im- 
poses constraints on the large and small x behavior of 
gp{x) and 5^(2;). Requiring that in the early /3- regime 
the e dependence should drop out of Xx(k,t), we find: 

Xx{k,t)^e l<t<r0. (16) 

Analogously, matching the a- and /3-regimes leads to: 

Xx{k, t) ^ e'^b-a),2a^b ^ ^ ^ (^7) 

interpolating between Xx ^ e^^/^ for t = Tp and Xx ^ 
e^^ for t = Ta, before decaying back to zero for t^Ta- 
All these results are valid close enough to the transi- 
tion but, as discussed previously, we expect crossovers 
as a functions of k and time. In order to study this is- 
sue numerically, we solve the full, wavevector-dependent 
mode-coupling equations for the self-intermediate func- 
tion Fs{k,t) for a dense colloidal suspension, which are 
directly coupled to the collective density fluctuations 
F{k,t) as 



dFs{k,t) 
dt 



Dok^F,ik,t) 



+ fdt' Ms{k, t - t')^^y^^ - 0, 

JO 



(18) 



dt' 



where Dq is the bare diffusion constant, and Ms{k,t) is 
the self-memory function that can be expressed as 



M,(fc,i) 



PqDq 
(2^) 



^ J dk' [k-k'c{k')y FA\k--k'\,t)F{k\t). 

(19) 

Here, po ~ N/V is the number density, k = k/|k|, c{k) 
is the direct correlation function, poc{k) = 1 — l/S{k). 
This equation must be solved simultaneously and self- 
consistently with mode-coupling equations for the full 
density fluctuations 

dnk,t),D.k^ f^^^, dF^^^^ 

Jo dt' 

(20) 



dt S{k) 
where 

PqDq 



M{k,t) 



jdk! |y(k,k')pF(|k-k'|,t)F(fc',t), 

(21) 




FIG. 6: Dynamic susceptibility x<p{k,t), Eq. (|22p . predicted 
by MCT for hard spheres at fixed volume fraction above the 
glass transition, (pc — (j) = 10~^ for various wavevectors from 
k — 19.35 to fc = 0.75 (in particle size units) from left to right. 
Power laws for the largest k are the asymptotic results and 
t'' with a = 0.312 and b = 0.583, while x<p{'t) ~ t describes 
well the data at small k (rightmost curve). 



and Vik,k') = k-k'c(fc') +k - (k-k')c(lk-k'l). These 
equations are solved for a model hard-sphere suspension 
with input from S{k) calculated from the Percus-Yevick 
closure at various volume fraction (jj. The wavevector 
cutoff is taken to be kc ~ 50 in units of the particle size, 
and the grid number is taken to be 100. The equations 
of motion are integrated with the algorithm of Fuchs et 
al. The induced susceptibility 



X4>{k,t) 



dFs{k,t) 



(22) 



is computed via numerical differentiation. In Fig. [5] the 
induced susceptibility is shown for different values of 
wavevector from those higher than the first peak in S{k) 
to those significantly below. The behavior of x<}>{k,t) 
for k close to the first diffraction peak displays the two 
power law regimes described in Sec. |lll When k is de- 
creased the power law describing how X4>ik, t) reaches its 
peak clearly shows an increasing value. This behavior is 
qualitatively compatible with the one discussed theoret- 
ically and found in numerical simulations of KCM's and 
atomistic liquids in (48j and experiments on granular ma- 
terials [l^l ■ It makes clear that corrections to the critical 
behavior are different depending on k. In particular, the 
limit fc ^ and T ^ Tc clearly do not commute. 

Note that, although the critical behavior is expected 
to be the same for X4>{k,t) and Xiit) in the NVE en- 
semble, the corrections to the critical scaling might not 
be the same. However, the same qualitative considera- 
tions regarding time and wavevector dependences should 
hold. Furthermore, using the bounds derived in I the 
present results yield direct quantitative predictions for 
the behavior of the dynamic susceptibility Xii^) in the 
hard sphere system since its scaling behavior is the same 
as to the one of x'^{k,t), which obviously follows from 
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Eqs. (nuiniiTaiTiD. 

III. DYNAMIC SUSCEPTIBILITY AND 
DIVERGING LENGTHSCALES IN KCMS 

A. Models and observables 

In this section we proceed with the computation of the 
dynamic susceptibihty xrit) and its comparison to the 
previously studied X4(0 kinetically constrained mod- 
els (KCMs) [l^. Our motivation here stems from the 
fact that the study of KCMs has greatly contributed to 
our present understanding of the dynamically heteroge- 
neous dynamics of supercooled liquids [H, Q, [H, 153, 
m, M, M, M, m, M,M, Ei. Moreover, the variety of 
available models allows one to grasp the variety of possi- 
ble behaviors that could possibly be encountered in real 
materials. Finally the relative simplicity of the models 
makes them suitable to large scale numerical simulations, 
which might help data analysis in real materials, while 
scaling laws and exact results can be obtained by stan- 
dard theoretical tools of statistical mechanics. 

KCMs are spin models (lattice gas versions also ex- 
ist [s^l) generically described by a simple, usually non- 
interacting Hamiltonian, and a set of dynamic rules with 
non-trivial constraints forbidding some of the transitions 
and therefore making the overall dynamics glassy. In the 
following we will focus on spin models characterized by 
the Fredrickson- Andersen Hamiltonian (5lj | , 

N 

H = Y,n^, (23) 

i=l 

where n,; = 0, 1 is a binary variable defined on each point 
of a hypercubic lattice in dimension d. Physically, = 
{ui = 1) represents a site i which is immobile (mobile), 
and has therefore an energy which is smaller (larger) than 
the average energy, given by 

= <T) = (24) 

The spins evolve with a single spin flip dynamics, so that 
the model dynamics is entirely defined by the transition 
rates between states 1 and 0, 

C. c 

= , = 1, (25) 

C. (l-c) 

where Ci is a kinetic constraint on site i which can become 
depending on the local environment of site i, therefore 
prohibiting some specific transitions. We shall study in 
detail two different spin facilitated models where the ki- 
netic constraint takes the following forms, 



and 

C, = 1 - - n,nk), (27) 

for 1- and 2-spin facilitated models, respectively. In the 
expressions for Ci the products are over nearest neighbors 
of site i. The constraints respectively become equal to 1 
when, respectively, at least 1 or 2 of their nearest neigh- 
bors is mobile, therefore capturing the idea of dynamic 
facilitation: mobile regions locally favor the creation of 
more mobility [l^, [5l| . 

Due to the presence of a heat bath, the dynamics of 
KCMs do not conserve energy. Physically this means 
that heat can be locally provided to a spin to allow the 
creation of a mobility excitation without the need to bor- 
row energy from the neighboring sites. In principle the 
results obtained from KCMs should then be compared 
to the NVT dynamics of molecular liquids. As opposed 
to the MCT results described above, no prediction can 
be made from KCMs concerning the role of a conserva- 
tion law for the energy. For kinetically constrained lattice 
gases, however, a quantitative comparison between spon- 
taneous fluctuations, X4(t), and fluctuations induced by 
a change of density, Xp(i), can be performed fisj . 

A second important consequence of the presence of a 
heat bath is that neither the fluctuation-dissipation rela- 
tion in Eq. ^ nor the inequality Eq. ^ apply to KCMs, 
and we are therefore left with three independent dynamic 
quantities, namely 

CpE{t)^N(6Pit)SeiO)), (28) 
X4{t) = N{5P\t)), 

where we have defined the instantaneous value of the 
energy density, 

1 ^ 

e{t)^j;jY.^.{t). (29) 

i=l 

Following earlier works on KCMs we choose to work with 
the persistence function as the relevant two-time dynam- 
ical object, 

1 ^ 

^W = ]^E^»W' (30) 

where Pi(t) denotes the persistence of the spin i between 
times and t. Its thermodynamic average, {P{t)), has 
recently been the subject of a number of theoretical stud- 
ies [H, 113, H m, [Hi, i^. Note that for KCMs, the 
Cauchy-Schwarz inequality could still be of some use. For 
the present variables this leads to 



C, = l-l[{l-n,), 

3 



(26) 



(31) 
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where cv{T) = dc/dT = e^/'^ /[T'^{1 + e^/^^f] is the spe- 
cific heat. The main difference between the inequalities 
(f3T|l and dH) is that the right hand side of ([3T|) is given by 
a correlation function which is not easily accessible in ex- 
periments, contrary to the susceptibility XT{i) appearing 
in (III). Of course, Cpsit) can be measured in numerical 
experiments, as shown below. 

At the level of the spatial correlations, two distinct 
correlators also need to be studied. 



ST{r,t) = (6P,{t)Sn,+r{Q)), 



(32) 
(33) 



Their Fourier transforms, S4{q, t) and Sxiq, t), can cquiv- 
alently be studied, and it is obvious that 5*4 (g = 0,t) = 
Xiit) and Sxiq = 0,t) — CpE{t), these quantities rep- 
resenting the volume integrals of the spatial correlations 
S4{r,t) and ST{r,t), respectively. 



B. Results for 1-spin facilitated FA models 

The one-spin facilitated FA model has been studied nu- 
merically and anal ytic ally in various spatial dimensions 
in much detail M, \M, H IH, [Hi, El, . These 

studies have shown that the model exhibits dynamic het- 
erogeneity and large spontaneous fluctuations of the two- 
time dynamics, although relaxation timescales grow only 
in an Arrhenius fashion as temperature is decreased. 



-A 



exp(A/r), 



(34) 



with A = 3 for d = 1 and A = 2 for d > 2. Interestingly, 
these works suggest that even strong material should dis- 
play dynamic heterogeneity. This was confirmed by sim- 
ulations and experiments [6^ which reported devia- 
tions from the Stokcs-Einstein relation, although the FA 
model itself presents no such deviations for d > 2. 

As usual, the four-point susceptibility Xii^) is found 
to have non-monotonic time dependence. Therefore it 
shows a peak, xtiT) = Xiit ^ Ta), whose position shifts 
to larger times and whose height increases when ternper- 
ature decreases. One finds dynamic scaling [isl. l47l. [ssj . 



X4 



exp(7/r) ^ rj/-^. 



(35) 



with 7 = 1 in all spatial dimensions. The correspond- 
ing spatial dynamic correlations have also been studied. 
Analytically, one can compute these quantities approxi- 
mately by making the assumption that the system can 
be described as an assembly of defects which diffuse in- 
dependently with diffusion constant D = c. This was 
called "independent defect approximation" in Ref. p^ . 
In three dimensions, one finds 



Siiq,t) ^ X4{t)S4q^e4it)], 
with a diffusively growing Icngthscale, 



(36) 



(37) 



and the scaling function 



S4x) = 2 



X- 1 



(38) 



Additionally the four-point dynamic susceptibility be- 
haves as follows. 



X4it) 



C2 l± 
2c \ T„ 



exp 



2t 

Try 



(39) 



with C2 a numerical factor. These predictions are in 
good agreement with direct simulations of the FA model, 
the only discrepancy being that the scaling function for 
S4{q,t) shows deviations from its 1/q^ predicted large q 
behavior when times become very large, t 3> Tq. 

The computation of XT{t) is easy given that the av- 
erage persistence obeys time temperature superposition, 
{P{t)) = f{t/Ta), the scaling function f{x) being well 
described, for times which arc not too long [60[, by 
a stretched exponential form, f{x) ~ exp{~x^), with 



(3 = 1/2 for d = 1 and f3 
immediately gets. 



1 for (i > 2. Therefore one 



/3 



A(3 / t 
XT(t) = -^j I — ) exp 



(40) 



This shows that XT(t) displays a non-monotonic time de- 
pendence with a peak arising at time t ^ Ta, diverging 

Xt ~ — 1/T^ when T goes to zero. Finally, the be- 
havior of Xt(^) before the peak, < ^ Tq,, is a power law, 
XT{t) ^ , P being the value of the stretching exponent 
characterizing also the a-relaxation. 

If one considers the quantity T^x^/cy appearing in 
the inequality (j4|); one finds at the peak. 



cv 



cxp(i/r) ^ xl 



(41) 



so that both sides of the inequality ^ have similar scal- 
ing properties at low temperatures in this model. This is 
not an obvious result given that these quantities arc not 
related by the thermodynamic relations and inequalities 
outlined in Sec. |T1 

Notice, however, that this similarity appears coinci- 
dental because the whole divergence of the first term in 
Eq. (|4T|) is due to the very strong temperature depen- 
dence of the specific heat at low temperature which itself 
results from the non-interacting FA Hamiltonian (j23p . In 
real materials, the specific heat is almost temperature in- 
dependent when the glass transition is approached and 
the growth of the term T'^xh'/cv is mainly due to the 
growing susceptibility XTit) itself. 

Following steps similar to those described in Ref. |T3| it 
is possible to compute both the correlator in Eq. ([5^ and 
its volume integral in Eq. (j28p within the independent 
defect approximation. In three dimensions, one finds for 
Sxiqjt) a scaling form very similar to Eq. (j36p . 



STiq,t)^CpEit)ST[q'm)] 



(42) 
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with 

Cp.«).c,(i-)oxp(-£li), (43, 

where Ci is a numerical factor, the corresponding corre- 
lation lengthscale 

Crit) = Vrt, (44) 
and the scaling function 

M^) - (45) 

These calculations show that, within 1-spin facilitated 
models, the physical content of the correlators 6*4 and 
St is essentially the same. Physically, this is because 
two sites are dynamically correlated, and therefore con- 
tribute to Si{r, t), if they are visited by the same diffusing 
mobility defect. Similarly, two sites contribute to Sxir, t) 
if one of them contains at time the first defect which 
will visit the second one for i > 0. This implies that the 
correlation lengthscales ^4 and both reflect the simple 
activated diffusion of point defects, and therefore con- 
tain the same physical information; Eqs. ((37)) and ([44]) 
show that they are indeed equal. Additionally, the spa- 
tial correlators S4{q,t) and Sriq^t) are found to differ 
in their detailed expression, Eqs. ([551 but they have 
the same asymptotic behaviors, ST{q^T ^ 1) '^-^ const 
and Sriq^T 3> 1) ~' 1/?^, reminiscent of an Ornstein- 
Zernike form. 

An additional piece of information derived from 
Eqs. (|40)) and ([43]) is the similar time dependence 
and scaling with temperature found for the quantities 
T'^XT{t) and CpE{t), despite the fact a fluctuation- 
dissipation relation such as Eq. ([3|) does not hold. Nu- 
merically we indeed find that both terms quantitatively 
differ, although merely by a numerical factor. 

The independent defect approximation is thought to 
be a good representation of the 1-spin facilitated model 
above its critical dimension, c? > 2, as confirmed by 
our numerical simulations in d = 3. In Fig. [71 we 
provide additional numerical evidence that these find- 
ings are also correct in d = 1 as well. The top figure 
shows that the time dependence and scaling with temper- 
ature of three different quantities, X4(t), Cp^{t)/ {T^cy) 
and T^Xt(0/cv are the same. Moreover, the Cauchy- 
Schwarz inequality ([3T|) is satisfied by our numerics, at 
it should be, while the inequality ^ derived for Newto- 
nian dynamics is found to be violated, by a factor which 
is about 2 at the peak. The bottom panel in Fig.[71shows 
Siiq, Ta) and Sxiq, Tq) for the FA model in d = 1. As pre- 
dicted from the independent defect approximation, both 
correlation functions are slightly different in their shape 
but share a common behavior, a plateau at small q and 
a 1/q^ decay at large q. We have checked that the equiv- 
alence ^4{t) ~ ^xit) also holds in numerical simulations, 
confirming that dynamic-dynamic and dynamic-energy 
spatial correlations are essentially equivalent quantities 




FIG. 7: Dynamic susceptibilities and spatial correlations in 
the one-spin facilitated FA model in one dimension. Top: 
Various dynamical susceptibilities are shown as a function 
of time for temperatures T = 1.0, 0.5 and 0.2 (from left to 
right). They behave similarly with time, and their peak scale 
shown with dots. Bottom: Dynamic struc- 
ture factors S4{q,t) and ST{q,t) at time t = Tc, for T = 0.3. 
Both functions behave as a constant at small g, and as 
at large wavevectors, as shown with a dashed line. 



in the context of non-cooperative KCMs. This is phys- 
ically expected since by definition of the kinetic con- 
straints in ([^S)) . it is those regions with high potential 
energy which trigger the dynamics of the nearby sites: 
this is the essence of the dynamic facilitation idea. 

We conclude this section on one-spin facilitated mod- 
els by briefly discussing the case of the East model [sot . 
which is defined with the same FA Hamiltonian ([23)) and 
is also a 1-spin facilitated model where the kinetic con- 
straint is defined similarly to Eq. ([^5)) . the only differ- 
ence being that the product appearing in (|26p is now 
restricted to only one neighbor in each spatial direction. 
This "hyper" -East model was called the North-or-East- 
or- Front (NEF) model in d 3 [H]. This direction- 
ality of the constraint makes the dynamics of the East 
model slower than that of the FA model, and relaxation 
timescales now grow in a super- Arrhenius fashion, so that 
the exponent A appearing in p4p becomes temperature 
dependent, A(T) ~ — In c(T). The East model is there- 
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fore a KCM for fragile glasses. Additionally, time tem- 
perature superposition does not hold. Relaxation is still 
described by stretched exponentials but the stretching 
exponent is also tem per ature dependent, with P{T) ^ T 
at low temperature [64l . [65j . Despite these qualitative 
differences between strong and fragile models, our main 
conclusions still hold. The three dynamic susceptibilities 
shown in Fig. [7] also track each other, and this is again 
the result of subtle compensations between the scaling 
of correlations functions and the strong temperature de- 
pendence of the defect concentration. Similarly, the two 
different lengthscales ^4 and also bear the same physi- 
cal content, although they now grow sub-diffusively with 
time [3] ■ This subdiffusive behavior affects the approach 
of the dynamic susceptibilities to their maximum. In 
the d — 1 East model, one finds that before the peak 
Xrit) t^'^'^K where the exponent b{T) « [3{T) should 
decrease slowly when T decreases. Wc find numerical 
values b sa 0.2 — 0.4 in the time window of our Monte 
Carlo simulations, where relaxation timescales increase 
from T„ - lO'^ to Ta ^ 10^. 



C. Results for a 2-spin facilitated FA model 

By comparison with one-spin facilitated models, much 
less is known about the behavior of 2-spin facilitated 
models, because relaxation does not proceed by activated 
diffusion (or even sub-diffusion) of point defects [111, . 
In some cases, asymptotic mechanisms have been de- 
scribed which show that relaxation timescales grow very 
rapidly when temperature is decreased, although no fi- 
nite temperature divergence is found [H, [s^ . In these 
mechanisms, relaxation occurs via the diffusion of "super- 
defects" whose concentration decreases when T decreases 
and whose size is itself an increasing function of tem- 
perature. For this reason these models are sometimes 
called "cooperative KCMs". Very recently, a KCM was 
specifically engineered to yield an example of a finite tem- 
perature singularity, but we do not discuss this example 
further [6^. 

An additional point of theoretical interest of coop- 
erative KCMs is that, when studied on Bethe lattices, 
they display a dynamical transition at finite tempera- 
ture which is reminiscent of the mode-coupling singular- 
ity described in Sec. [Ill Moreover, dynamic fluctuations 
can be studied in some analytic detail in the Bethe limit, 
while no analytic study of dynamic fluctuations on finite 
dimensional lattices is available. 

We now focus on the 2-spin facilitated model in di- 
mension d ~ 2, the "22FA model", as a specific exam- 
ple of a cooperative model. Our choice is motivated by 
the relatively large number of earlier studies dedicated to 
this model [5l|, [S: 
also considered 



5^ . the fact that its Bethe limit was 



and that it is sufficiently far from 
its mean-field limit that deviations from mean-field be- 
havior are clearly observed. It was indeed realized early 
on that the model docs not display a power law diver- 



gence of its relaxation time at finite T [5ll . |52| | , contrary 
to more constrained models where numerics seemed to 
indicate the presence of a mean-field like singularity [H^] , 
now discarded [H, [IJ ■ 

Adapting the general results of Ref. [s^ to the spe- 
cific example of the 22FA model, we expect the following 
scaling results. The relaxation time grows as 



exp 



(?) 



exp 



a exp [ ^ 



(46) 



where a is a numerical factor. The double exponential 
divergence makes the 22FA a very fragile glass-former 
model. 

The scaling of the four-point dynamic susceptibility is 
obtained as follows. At a given temperature, relaxation 
occurs via the diffusion of super-defects of size i{T). By 
coarse-graining the system up to size £, relaxation then 
resembles the diffusion of independent defects, and the 
results of the independent defect approximation can be 
carried out. Therefore, we expect xt ^ cj^, where ce is 
the concentration of super-defects. Using the results of 
Ref. [5^, we get 



X4 ~ exp 



(!) 



(47) 



We evaluate the leading divergence of xt(0 by as- 
suming time temperature superposition, i.e. XT{t) = 
df{t/Ta)/dT. Using ^ we get Xt ^ exp{a/c)/{T^c), 
up to an irrelevant numerical prefactor. As a conse- 
quence, the right hand side of the inequality ([4]) scales 
as 



rp2 

Cv 



(IriTa) 



(48) 



By comparing Eqs. p5)) and ([47]), we conclude that 
the dynamic heterogeneity quantified through X4(i) and 
r^XT(*)/cy are very different, since Xiit) is predicted to 
diverge as a power of r^, while the term involving xt 
should diverge only logarithmically with r^. For coop- 
erative models, the "coincidental" compensation due to 
the specific heat arising in non-cooperative model is not 
effective. 

Since these results are expected to hold only very close 
to T = 0, we have performed numerical simulations of the 
22FA model. In these Monte Carlo simulations, we cover 
the temperature regime T = 2.6 down to T = 0.43, which 
corresponds to about 7 decades of relaxation timescales. 
In this temperature window, Tq, cannot be fitted with an 
inverse power law Tq, ~ (T — Tc)~" as in the Bethe limit, 
showing that strong non-mean-ficld effects are indeed 
present. However, the form (j46p is not completely suc- 
cessful either, suggesting that the true asymptotic regime 
is beyond the realm of numerical simulations (see [6^ for 
a discussion of this point in a similar context), and that 
the numerical regime lies somewhat in a crossover regime. 

In Fig.[5]wc compare the evolution of the peak of X4(t) 
and the corresponding peak in T'^XtW/'^v ^^^r the entire 
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FIG. 8: Dynamic susceptibilities and spatial correlations in 
the 2-spin facilitated FA model in two dimensions. Top: Com- 
parison between the peak values of X4 and T^Xt/cv for dif- 
ferent temperatures covering about 7 decades in relaxation 
timescales. The full line represents the proportionality be- 
tween both quantities. Bottom: Dynamic structure factors 
S4{q,t) and S'T((7,i) at time f = ~ lO'^ for T = 0.428. 
Both functions behave as a constant at small q, but have dif- 
ferent large q behaviors since the 1/q^ dashed line is consistent 
with S4, only. 

temperature range we have been able to access numeri- 
cally. Quite strikingly wc find that both functions scale 
very similarly on the whole temperature range. A sim- 
ilar result was recently reported for a cooperative con- 
strained lattice gas in two-dimensions [i^. This similar- 
ity holds also at the level of the whole time dependence 
(not shown). From numerical simulations only, we would 
therefore conclude that the coincidence between the two 
terms already found for non-cooperative models also ap- 
plies in cooperative models. This numerical evidence is 
contradicted by the asymptotic analytic arguments given 
above. A possible way to reconcile these results is to as- 
sume that the temperature regime we have studied in 
the simulations is still too close to the mean-field Bethe 
lattice limit, where the scaling X4 ~ T'^Xt/cv is indeed 
expected to hold. This argument is however clearly weak- 
ened by the fact that many obscrvables (timescales, per- 
sistence functions, and others, see Ref. [60]) show visible 
deviations from their mean-field limit in the same tem- 



perature regime. 

In the bottom panel of Fig. [51 we also show the com- 
parison of the spatial correlations (|32p and ([55)) measured 
in Fourier space. Whereas both correlators were found 
to be very similar for non-cooperative models, numerics 
clearly reveals that the shapes of the dynamic structure 
factors 5*4 and St differ. While 5*4 ~ l/q^ seems to hold 
at large wavevectors, 9^4 S> 1, we find a different behav- 
ior for St, namely St ~ \/q^'^. Note that this fit is not 
very satisfactory, revealing a more complex structure of 
this correlator, possibly related to the presence of two 
Icngthscales in the model: the size of the super-defects, 
and the typical distance separating them. We conclude 
that dynamic-dynamic and dynamic-energy correlations 
might contain slightly distinct physical information in co- 
operative KCMs. This is physically expected because an 
isolated defect, which represents a positive local fluctua- 
tion of the energy, cannot diffuse and relax the neighbor- 
ing sites. Therefore the correspondence between energy 
fluctuations and dynamical fluctuations is not one-to-one 
as in, e.g., the 1-spin FA model. By this argument one 
can predict that >S't(?', t) < Si{r,t) at small r, and there- 
fore a faster initial decay of STir^t) with r. In Fourier 
space, this means a slower large q decay of ST{q,t) than 
that of S4{q, t), as observed in Fig. [5| 

Nevertheless, a dynamic correlation lengthscale can be 
defined from both 5*4(5, i) and STiq,t) as the inverse of 
the wavevector above which structure factors start to de- 
cay. The data shown in Fig. [5) clearly indicate that these 
two lengthscales are very close. A possible interpretation 
is that despite their complex structure, super-defects re- 
main associated with some positive energy fiuctuations, 
so that the lengthscale extracted from three- and four- 
point functions could indeed be equivalent, as in the case 
of non-cooperative models. A similar situation was en- 
countered in our atomistic simulations in I. 



D. Remarks and open questions on ensemble and 
dynamics dependence and KCMs 

We have studied in the context of kinetically con- 
strained spin models (KCMs) the dynamic susceptibil- 
ity XTit) and the associated three-point dynamics-energy 
spatial correlations, and their link with the more stan- 
dard four-point susceptibility Xii^)- Although the ther- 
modynamic relations derived in I for supercooled liquids 
do not hold for kinetically constrained spin models (be- 
cause energy is not dynamically conserved), they seem to 
be approximately valid. 

The underlying reason is that in non-cooperative 
KCMs the energy fiuctuations that are important for the 
dynamics are effectively conserved because of the kinetic 
constraint. This is clearer in the example of the 1-spin fa- 
cilitated FA model where a facilitating spin can disappear 
only by annihilation with another facilitating spin. Simi- 
larly, a facilitating spin can be created only by branching 
from another facilitating spin. But these two processes 



14 



happen very rarely (see Refs. [13, [H, for a detailed 
analysis and discussion of timescales). Therefore, the 
main relaxation mechanism is diffusion of the facilitat- 
ing regions (energy fluctuations) which are conserved in 
an effective way, as assumed in the independent defect 
approximation. 

Comparing our results for KCMs to the general the- 
oretical considerations of I opens interesting issues re- 
lated to the applicability of KCMs to supercooled liq- 
uids. Since dynamical fluctuations strongly depend on 
statistical ensembles and microscopic dynamics, this im- 
mediately raises important questions: 

• For which ensemble arc the dynamical fluctuations 
of real liquids supposed to be described by KCMs? 

• What type of liquid dynamics should one choose to 
compare real dynamical fluctuations to the predic- 
tion of KCMs? 

These questions are clearly related to the coarse- 
graining procedure that is often invoked [5l|, [H, [tOI, 
but never truly performed, to map real liquids to KCMs. 
Were this procedure known, the answer to the previous 
questions would be clear. Unfortunately, this formidable 
task has not yet been accomplished. On the other hand, 
our results show that this issue is important if one wants 
to compare KCM predictions for dynamic susceptibili- 
ties Xi s-nd XT to experimental and numerical results on 
realistic models. 

For kinetically constrained spin models, the answer to 
the first of the above questions seems fairly easy even 
without the coarse graining procedure. Only in the most 
general ensemble where all conserved quantities fluctu- 
ate does one have lim^^o 'S'4 (9, i) = 5*4(0, t). Since this 
equality holds in KCMs, we conclude that KCMs should 
apply to real liquids in the most general statistical en- 
semble, i.e. NPT for most practical purposes. 

The second question is instead much more subtle. 
From a general point of view since there is no conserved 
quantities in spin models, KCMs could be thought as rep- 
resentative of a dynamics without conserved quantities. 
Of course all physical dynamics should at least conserve 
density. However, if one considers Brownian dynamics 
for supercooled liquids for which temperature is the rele- 
vant control parameter, while density plays a minor role 
(see section II. E. 3 in I), it might be reasonable to expect 
that density fluctuations do not couple strongly to dy- 
namical fluctuations. One is then tempted to conclude 
that KCMs are models of real liquids with Brownian or 
stochastic dynamics. 

However, this tentative answer is contradicted by sev- 
eral facts. First, real supercooled liquids obviously evolve 
with Newtonian dynamics. Second, we just discovered 
that the inequality (|3ip provides a good approximation 
to Xiit) for KCMs. A similar result holds for liquids with 
Newtonian dynamics in the NPT ensemble (see I) but 
not for liquids with stochastic dynamics (Tlj . 

Taking the opposite view that KCMs represent, for 
some unclear reason, liquids with Newtonian dynamics 



is also unsatisfactory because the saturation of the in- 
equality pTjl in KCMs is principally due to the behavior 
of the specific heat that decreases exponentially fast as 
temperature decreases. But a very small specific heat 
is incompatible with experimental measurements of the 
thermodynamics of supercooled liquids [ill . Correcting 
for this fact as in Ref. [zl] then leads to poor estimates 
of X4(0 dynamic response functions, in disagreement 
with atomistic simulations p^ . [l9| . 

The case of kinetically constrained lattice gases is less 
problematic if taken as models of glass/jamming transi- 
tion in hard sphere systems, rather than molecular liq- 
uids. In this case, the only conserved quantity that mat- 
ters is the density and therefore there are no ambiguities 
since density is conserved both in kinetic lattice gases 
and in real systems. 

KCMs provide a natural mechanism explaining corre- 
lations between energy fluctuations and dynamic hetero- 
geneity. However, in order to compare even qualitative 
predictions of KCMs with experimental or numerical re- 
sults for dynamical fluctuations, one has to understand 
clearly in what ensemble and for what dynamics KCMs 
predictions hold. This certainly highlights the impor- 
tance of a microscopic derivation and more detailed jus- 
tification of KCMs. 



IV. NUMERICAL RESULTS FOR TWO 
MOLECULAR GLASS-FORMERS 

A. Models and technical details 

In this section we report our numerical calculations of 
the dynamic susceptibility xt(0 in two molecular glass- 
formers which have been extensively studied in numerical 
simulations: a binary Lennard-Jones (LJ) mixture [t^ . 
considered as a simple model system for fragile super- 
cooled liquids 0, and the Beest, Kramer, and van San- 
ten (BKS) model, which is a simple description of the 
strong glass- former silica [giI . pTSj . For both models we 
have investigated the behavior of the dynamical fluctua- 
tions performing microcanonical simulations at constant 
energy, E, number of particles, iV, and volume, V , by 
solving Newton's equations of motion [t^. For the LJ 
system we have also simulated two types of stochastic 
dynamics, namely Brownian and Monte-Carlo dynam- 
ics [zl. 

We follow the dynamical behavior of the molecular liq- 
uids through the self-intermediate scattering function, 

j.^(k,i)^^_L|^e^k.[r,(*)-r,(o)]^^ (49) 

where the sum in Eq. runs over one of the species 
of the considered hquid [A or B in the LJ, Si or O 
for silica). We denote by /s(k, t) the real part of the 
instantaneous value of this quantity, so that we have 
F,(k,t) = (/,(k,i)). 
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As usual, the four-point susceptibility. Xi{t)j quantifies 
the strength of the spontaneous fluctuations around the 
average dynamics by their variance. 



Xi{t)^N^[{f!iKt))-F^{k,t)] 



(50) 



In principle, Xii^) Eq- (|50p retains a dependence on 
the scattering vector k. Since the system is isotropic, we 
circularly average (HI]) and ([5D|) over wavevectors of fixed 
modulus, and we mainly consider results for |k| = 7.21 for 
the LJ system, and |k| = 1.7 A"^ for the BKS. These val- 
ues respectively represent the typical distance between A 
particles, and the size of the Si04 tetrahedra. Finally, we 
use finite difference to evaluate the temperature deriva- 
tives involved in 



XT(i) = ^F,(k,t). 



(51) 



We have given an extensive account of the models, nu- 
merical details and parameters used in 1. Therefore, we 
refer readers interested in the technical details concerning 
the simulations to 1 [ll]. Note also that we will neglect, 
as it is justified in I, the role density fluctuations on dy- 
namical correlations. Therefore we will focus on Xi^"^ 
(instead of Xi^'^) £^ud xt obtained deriving with respect 
to temperature at fixed volume (and not fixed pressure). 



B. Time dependence of dynamic susceptibilities 

1. Time behavior of xrit) 

Our results for the dynamic susceptibilities, XT{t), are 
presented in Fig. [H] for both the LJ and BKS models. 
For a given temperature, the qualitative time depen- 
dence of xt(0 observed in Fig. [S] resembles the one al- 
ready reported for X4(i): Xt(0 presents a peak for a 
timescale close to Tq. This is very natural since by def- 
inition XT{t = 0) = XT{t — *■ oo) = 0, and it is for times 
i « Tq, that the dynamics is most sensitive to temper- 
ature changes. We have shown the quantity |xt(OI 
these figures, as XT{i) is obviously a negative quantity: 
raising the temperature makes the dynamics faster, and 
hence two-time correlators smaller, so that dFg/dT < 0. 

More quantitatively, we expect the two-timescale re- 
laxation of the averaged dynamics to lead to a com- 
plex time behavior of xt(Oi similar to that predicted 
for Xi{t) [III. Within MCT, we expect (see Sec. Ill two 
distinct power laws, xt ~ t"^ followed by xt ~ t^, to 
describe the approach to the maximum of xt, the ex- 
ponents a and b being already constrained to the val- 
ues they take when fitting the averaged dynamics using 
MCT. From the study of KCMs only the approach to 
the peak can be predicted since the short-time dynamics 
contains no clear relaxation towards a plateau due to the 
coarse-grained nature of the models [H, HI, [6l|. Again, 
a power law approach to the peak is expected. 

In Fig. [S] we compare our numerical results for XT(t) 
to power law behaviors shown as dashed lines. On the 




FIG. 9: Normalized xt as a function of time for various tem- 
peratures in a binary Lennard- Jones mixture (top) and the 
BKS model for silica (bottom), obtained from molecular dy- 
namics numerical simulations. LJ: T = 2.0, 1.0, 0.75, 0.6, 
0.5, 0.465, 0.432 from left to right. BKS T = 6000, 4650, 
4000, 3550, 3200, 3000, and 2730 K from left to right. We 
have taken the absolute value since xt is a negative quantity. 
Power law fits of the time dependence are discussed in detail 
in Sec. lIVI The value of the exponents at short and long times 
are 0.32 and 0.45 in (c), 0.3 and 0.5 in (d). 



restricted time window of the simulations there is obvi- 
ously some freedom in the fitting procedure so the ex- 
ponents we report should be considered as an empiri- 
cal quantitative description of the true time dependence 
of these functions. As discussed already in the case of 
X4(i) [111, corrections to the asymptotic scaling laws de- 
rived by theoretical approaches should be expected in the 
reduced time regime of the molecular simulations. In the 
LJ system we find that the time behavior of XT(i) can be 
described by the exponents a w 0.32 and b w 0.45 with 
the tendency that these exponents very slowly decrease 
when T decreases. For the BKS system we find a simi- 
lar quality of the fits with a ~ 0.3 and b w 0.5 with no 
systematic dependence in temperature. 

The values of these exponents compare reasonably well 
with the MCT predictions obtained above. For the LJ 
system, the von-Schweidler exponent is estimated to be 
b « 0.51 from fitting the averaged dynamics in the /3- 
relaxation regime |74l |. while direct computations pre- 
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diet b = 0.62 [t^]- Both values are close to our finding, 
b « 0.45, although they both shghtly overestimate it. 
The exponent a describing the dynamics in the early /3- 
regime was not directly fitted, but using the known rela- 
tions between MCT exponents its value is predicted to be 
a = 0.29 (for b = 0.51) and a = 0.32 (for b = 0.62). This 
is again consistent with our finding for xt(Oi ~ 0.32, 
in this time regime. From the point of view of MCT, we 
suggest that focusing on xt is a more powerful way to di- 
rectly measure the exponent a (this might be interesting 
from an experimental point of view as well). Finally for 
BKS, fitting of the average dynamics provides the value 
b = 0.62, from which a = 0.32 is deduced from known 
MCT relations [HH. These two values again compare 
relatively well with the time behavior found for XTit), 
namely a ~ 0.3 and b w 0.5. 

Applying results from KCMs to real liquids, one would 
predict the time dependence of XT{t) when approaching 
the peak to be XT{t) ~ t for an Arrhenius liquid mod- 
elled by the 1-spin facilitated model in three dimensions, 
while XT{t) ~ t''^^-' is predicted for fragile liquids mod- 
elled by the East model. Our numerical results for BKS 
silica are not consistent with the FA model predictions 
and are, quite unexpectedly, more compatible with the 
smaller exponents observed in the fragile East model re- 
ported in Section UlI Bl The small b{T) exponents of the 
East model compare however well with the behavior of 
Xxit) found in the LJ system. In particular, the fact that 
b{T) decreases with decreasing T is correctly predicted by 
fragile KCMs, as opposed to the constant b predicted by 
MCT. For a summary of these results, see Table I. 



2. Comparison between X4(^) ci^d XT{t) 

It is interesting to compare the exponents found nu- 
merically for XT{t) to the ones of Xi{t) measured in 
the NVE ensemble for Newtonian dynamics since the- 
ory predicts some relations between them. The latter 
exponents were already studied in Ref. [l3| for the LJ. 
Numerically no power law behavior X4(0 ^ is found 
in the short-time behavior of X4(t) in the Newtonian dy- 
namics of both the LJ and BKS systems. This is due 
to the fact that thermal vibrations strongly affect the 
short-time dynamics of these liquids. Two power-law 
regimes are however clearly observed in the stochastic 
simulations where phonons are either overdamped (Brow- 
nian dynamics), or absent (Monte-Carlo dynamics). Our 
Monte-Carlo results for Xi{t) i^i the LJ are presented in 
Fig. [TO] (top) where we have fitted the early and late (3 
regimes with two power laws with exponents a w 0.37 and 
b « 0.7, respectively. For the BKS we performed New- 
tonian dynamics simulations only. Hence, we only have 
results on the exponent b from Xi measurements, which 
is found to increase from 0.65 to 0.85 upon lowering the 
temperature: this is an opposite behavior compared to 
the LJ where b decreases. This might suggest a different 
temperature behavior of b in strong and fragile liquids. 
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FIG. 10: Top: Four-point susceptibility x^i^) in the binary LJ 
mixture with Monte-Carlo dynamics for T — 2.0, 1.0, 0.75, 
0.6, 0.5, 0.47, 0.45 and 0.43 (from left to right), the lowest 
temperature being highlighted with open circles. Power laws 
X4 ~ and Xi ~ are indicated with dashed lines in 
the early and late /3-regimes, respectively. Bottom: X^i^) 
is shown for T = 0.45 for NVE Newtonian, Brownian and 
Monte Carlo dynamics as a function of a rescaled time chosen 
so that all Xi'^ overlap near the alpha relation. We chose 
t = t for NVE Newtonian dynamics, t — f/24 for Brownian 
dynamics, t = t/100 for Monte Carlo dynamics. No rescaling 
of the vertical axis is performed: The agreement between the 
3 types of dynamics is remarkable. 



This trend is partly captured by KCMs. 

MCT predicts that xrit) and Xi^^{t) have the same 
critical scaling. KCMs predictions are ambiguous so we 
follow the numerical results obtained in Sec. IIIIl i.e. 
X4(0 ~ Xt(^)- I'^ both LJ and BKS systems, the expo- 
nent a is the same for both susceptibilities, as predicted 
by MCT. The results for b are more difficult to interpret: 
although b for Xi is systematically larger than for XT^ 
the ratio between the two exponents is not 2 either, so 
that neither MCT nor KCMs approaches really describe 
this aspect of our numerical results. For a summary of 
these results, see Table I. 

What comes nicely out of the simulations, however, 
is the fact, predicted on general grounds in I and 
within MCT above, that NVE Newtonian, Brownian 
and Monte Carlo dynamics display similar time depen- 
dences for the dynamic susceptibility Xi (t) ■ This is strik- 
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Observable 


LJ 


BKS 


MOT (LJ) 


MOT (BKS) 


KCM (IFA) 


KCM (East) 


e (xt) 


0.33 


0.14 


0.43 


0.43 


0.25 


oc T 




0.39 


0.18 


0.43 


0.43 


0.5* 


oc T 


a (xt) 


0.32 


0.3 


0.29-0.32 


0.32 


NA 


NA 


a ixr^) 


0.37^ 


NA 


0.29-0.32 


0.32 


NA 


NA 


b (xt) 


0.45 


0.5 


0.51-0.62 


0.62 


1 


P{T) oc T 


b ixr"^) 


0.7 


0.65-0.85 


0.51-0.62 


0.62 


2* 


P{T) oc T 



TABLE I: Summary of the different results for exponents 0, a and b, describing the peak amplitude and the time dependence 
of T\xt\/ y'^cy and x^^^ (see text). NA: not applicable; ^: Obtained from MC dynamics; *: Ambiguous - do KCMs describe 
X4^^ Newtonian or Xi^'^ Brownian (= '^'')? 



ingly illustrated in Fig. [10] (bottom) which shows Xiit) a-t 
a single temperature, T — 0.45. The results for the three 
dynamics almost perfectly overlap for timescales larger 
than the plateau regime in Fg (k, t) . 



C. Peak amplitude of dynamic fluctuations 

We now focus on the amplitude of the peak observed 
in the various susceptibilities. In Fig. 1111 we present 
our numerical results for Xa^^i T'^XtI^-v a-nd their sum 
X^^^ obtained from the Newtonian dynamics of both 
the LJ and BKS models. When temperature decreases, 
all peaks shift to larger times and track the a-relaxation. 
Simultaneously, their height increases, revealing increas- 
ingly larger dynamic fluctuations as the glass transition 
is approached. 

The main observation from the data displayed in 
Fig. [Til already made in Ref. and in I, is that in 
both LJ and BKS systems the term T^x^/cv while be- 
ing small; ^ O(10~^), above the onset temperature of 
slow dynamics, grows much faster than Xi^^ when the 
glassy regime is entered. As a consequence, there exists 
a temperature below which the temperature derivative 
contribution to the four-point susceptibility Xi^'^ dom- 
inates over that of Xa^^ ■ This crossover is located at 
T « 0.45 in the LJ system, T « 4500 K for BKS sil- 
ica. Remarkably, the conclusion that T^Xr/cy becomes 
larger than Xi^^ at low temperatures holds for both 
strong and fragile glass- formers. Experimental and the- 
oretical consequences of this observation were discussed 
in Refs. 

Following Ref. [i^l we have chosen to present the evo- 
lution of the amplitude of the dynamic susceptibilities 
as a function of Tq rather than T because it is in this 
representation that dynamic scaling might emerge. For 
the LJ system we find that all susceptibilities can be de- 
scribed by power laws, % ^ r^, in some intermediate, and 
therefore subjectively defined, temperature regime with 
following exponents: k, 0.39 for Xa^^ ^ ^ ~ 0.46 for 
and 9 w 0.67 for T'^xl/cv- For the BKS model, 
we find 9 ~ 0.27 for both Xa^'^ and, in a more restricted 
time window, T^Xt/cv while we find ~ 0.18 for Xa^^ ■ 

The theoretical considerations given above show that 




Ta (ps) 

FIG. 11: Various susceptibilities in the binary LJ mixture 
obtained from the A particles dynamics (top) and the BKS 
model for silica from the Si ions dynamics (bottom). Dashed 
lines indicate power law behavior with exponents 0.46, 0.39 
and 0.67 (from top to bottom in the LJ system), and 0.27 and 
0.18 (from top to bottom in the BKS model). In all cases, 
T^Xt/cv is smaller than Xa^^ high temperature, but in- 
creases faster and becomes eventually the dominant contribu- 
tion to Xa^"^ in the relevant low temperature glassy regime. 



these exponents should be related, within MCT, to the 
exponent 7 describing the divergence of Tq, close to 
Tc. The prediction is that 6* = I/7 for xf^^ and 
Txt/^ m while 9 = 2/7 for Xa^^ and T^xl/cy 
Fitting of the relaxation times has shown that 7 w 2.35 
for both LJ and BKS systems, so the exponents 1/7 = 
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0.426 and 2/7 = 0.851 should be observed in Fig. [TTl 
The exponent for x^^^ is reasonably well described by 
MCT predictions in the LJ system, an agreement already 
reported in Refs. [s^, [zH (see Table I). The agreement 
deteriorates somewhat for Txt I \fcv ■ The MCT predic- 
tions fail however strongly in the BKS system, for which 
the value 0.18 is found instead of the expected 0.426 for 
y^^^ and xt, although in a temperature regime where 
Arrhenius behavior is already observed. No clear power 
law can be seen in the mode-coupling regime seen in [6l| . 
In principle, the behavior of xt(0 is completely tied to 
the one of the average two-time correlators already stud- 
ied in [6l|, but xt(0 provides a more detailed analy- 
sis of the dynamics with no fitting procedure required. 
Therefore the failure of MCT to capture the behavior of 
Xt(0 suggests that MCT, despite the claims of [6l|, docs 
not satisfactorily describe the dynamical behavior of this 
strong glass-former. 

Finally we find that T'^Xj^/cy and Xa^'^ behave some- 
what differently in the temperature regime where power 
law fits are performed. This is not surprising. We have 
extensively discussed in I the fact that simulations are 
typically performed in the relatively high temperature 
regime where both terms contributing to Xi^"^ a-re com- 
parable. Since they are predicted to have different scal- 
ing behaviors, the intermediate value for the exponent 6 
reported for Xi^'^ simply results from this crossover. 

The power law regimes we have discussed do not de- 
scribe the whole temperature range studied for the LJ 
system. For T < 0.47 the growth of all dynamic sus- 
ceptibilities with Tq, becomes much slower, perhaps loga- 
rithmically slow, but we do not have a sufficient range of 
timescales in this low temperature regime to draw more 
quantitative conclusions. We have moreover checked that 
this saturation is not the finite size effect expected if fiuc- 
tuations are computed in too small a system size (67j . 
see I. Interestingly, no such saturation can be observed 
in the BKS system. Therefore we do not know how to ex- 
trapolate the present numerical results towards the glass 
transition temperature, and compare our simulations to 
the result, reported in Ref. [l^ . that dynamic susceptibil- 
ities have typically the same value at Tg for liquids with 
very different fragilities. We can simply state from our re- 
sults that this fragility-independence cannot hold at all 
temperatures since Fig. [TT] clearly shows that dynamic 
susceptibilities grow at different rates in different sys- 
tems. We are currently investigating this point in more 
detail 

The saturation of the LJ dynamic susceptibilities ob- 
served at low T seems consistent with the theoretical 
expectation [H, [H, [g^jT^UMJS^IllI , and the exper- 



imental confirmation 1^, [S^, [SJ, 13 ^^at dynamic fluc- 



tuations and lengthscales grow very slowly when T is de- 
creased towards Tg. From the fragile KCMs perspective, 

one would for instance expect that X4 ~ Ta^^'^ with an 
exponent 6{T) which decreases linearly with T [63l. 
while logarithmic growth, Y4 _~ (log 70)^1 is predicted by 
activation based theories [3i 



CONCLUSIONS AND PERSPECTIVES 



This paper describes the second part of our investiga- 
tions of dynamical susceptibilities started in I [T6j . In this 
second work we have illustrated the general conclusions of 
I by making explicit the predictions of MCT and KCMs 
concerning spontaneous dynamical fluctuations (encoded 
in Xiii)) and induced one (given by Xx{i))- These theo- 
ries predict the detailed dependence of these two quan- 
tities both as a function of time and of temperature (or 
density). As discussed in I, special care must be devoted 
to the choice of statistical ensemble and microscopic dy- 
namics, with the rather spectacular prediction of MCT 
that X4(^) should coincide (or at least display the same 
scaling) for Newtonian dynamics in the NVE ensem- 
ble and for Brownian dynamics in the NVT ensemble, 
but differ from the result for Newtonian dynamics in the 
NVT ensemble. The predictions coming from KCMs are 
much less clear about this particular point, since there 
is some intrinsic ambiguity about which ensemble and 
which dynamics these models are supposed to describe. 

We have compared these predictions with numerical 
simulations of models of supercooled liquids. Overall, 
as shown in Table I, MCT fares reasonably well at ac- 
counting for the detailed shape of Xi{i) and XT{t) of the 
Lennard-Jones system, in a restricted temperature re- 
gion where MCT can be applied. As for the values of 
the exponents, our aim was to present a rather qualita- 
tive comparison focusing more on compatibility than on 
precise tests, which arc beyond the scope of this work, 
and probably of MCT as well. Quite remarkably, the 
exponents used to fit these higher order correlations are 
indeed compatible with those measured on two-point cor- 
relation functions, with quantitative variations that can 
perhaps be attributed to preasymptotic effects. Further- 
more, the predicted ensemble dependence of these quan- 
tities is very clearly highlighted by our numerical results. 
We have also shown that the wave-vector dependence of 
Xi{i) can be quahtatively accounted for within MCT. On 
the other hand, the features of the dynamical suscepti- 
bility of the BKS model for the strong silica glass are not 
quantitatively well explained by MCT. Similarly KCMs 
fail to describe quantitatively the results obtained in the 
BKS model, but the systematic temperature dependence 
of the exponents describing Xiit) appears somewhat nat- 
ural from this perspective. 

Among open problems, we should primarily emphasize 
the major problem of extending MCT to allow for acti- 
vated events. A detailed prediction of X4(0 and of the ge- 
ometry and exponents of dynamically correlated regions 
in the deeply supercooled region would be important to 
compare with future experiments (see [13, Hi], [s^ for pre- 
liminary elements in that direction). The generalization 
of these predictions to the aging regime would also cer- 
tainly be relevant to analyze the cooperative dynamics of 
deeply quenched glasses. 
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Appendix A. DYNAMIC SUSCEPTIBILITIES IN 
THE p-SPIN MODEL 

A.l. General discussion and results 

Much intuition concerning dynamic heterogeneity has 
been gleaned from the study of mean-field spin-glasses. 
In particular, Franz and Parisi first pointed out that a 
quantity analogous to X4(^)j which can be computed ex- 
actly in mean- field p-spin models, should show non-trivial 
features [llj , which prompted the study of dynamic fluc- 
tuations in simulations of atomistic glass-forming liq- 
uids [l^ . The growth of a dynamic susceptibility in this 
model was properly interpreted in terms of a growing dy- 
namical length scale, which diverges at Tc- The same sce- 
nario, complete with a temporal behavior of X4(0 iden- 
tical to that in the p-spin models, exists in mean-field 
models that have no underlying thermodynamic critical 
point [H, [S^ . It should also be noted that this scenario 
is perhaps more general than appreciated, since it ap- 
pears to also exist in models on compact lattices with no 
quenched disorder and short-ranged interactions, at least 
in the limit of large dimensionality [8 71 . an d models with 
long-ranged, Kac-like interactions [8^, [S^] ■ 

Applying the above diagrammatic analysis to p-spin 
models for which no conserved quantities exist, one finds, 
in agreement with BB, that X4(i) is determined by lad- 
der diagrams only. Hence, its critical behavior has to 
be the same as the dynamical response XTit) and is 
given by Eq. (fTTj) . Similarly the susceptiblity xfp(^) 
introduced by Franz and Parisi is found to follow the 
same scaling behaviour. As discussed below, Franz and 
Parisi [ll[ study the quantity XFp{t) = dC{t)/de, where 
C{t) = {si{t)si{Q)) and e is an infinitesimal field cou- 
pling the system's configuration at time t to its initial 
state at time 0. Using linear response theory they ar- 
gue that dC{t)/de and X4(0 ^-^e equal. We find in- 
stead that dC{t)/de is equal to the sum of x^i^) ^-nd an- 
other non- vanishing contribution. However dC(t)/de — 
XFpit) = N-'J:^J^ dt'{s,{t)smsjit')s,{t'+)), where 
Si{t) are the response field. Hence, it is given by ladder 
diagrams similar to the ones contributing to X4(0- Thus 
we expect that XFp{t) and Xiit) behave similarly close 
to the critical point. 

In the following, we shall present a careful numerical 
comparison between the dynamic susceptibility XFp{t) 
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FIG. 12: Time dependence of the dynamic susceptibilities 
XT(t) (thick lines) and Xfp(^) (thin lines) in the p = 3 mean- 
field p-spin model for temperatures approaching Tc from 
above. Note the wide range of timescales covered in this 
graph. From left to right, {T -T^)/T^ = 10"^ 10"'', 10"^ 
10~*. The asymptotic power law regimes are shown as dashed 
lines. The susceptibilities grow as , and in the micro- 
scopic, early and late beta regimes, while the height of the 
maxima scale with their as x* ~ t^^'^ . For p = 3, one has 
a = 0.395, fe = 1 and 7 = 1.765. 

and XT{t) integrating the intcgro-diffcrcntial equations 
derived in (l]| for p-spins models. This comparison de- 
cisively confirms the previous analytical results. A much 
smaller time window was studied in (33 | , and it was not 
clear that asymptotic regimes had been observed. 

One technical difficulty is that it is numerically diffi- 
cult to calculate xfp(0 very close to Tc- Here, we modify 
the method developed by Kim and Latz for the aging p- 
spin model [93| to accurately integrate the equations on 
Xfp(0 derived in Ref. [ll| much closer to Tc than has 
been reported in previous work. The dynamical equa- 



tions arc presented in [Appendix A A. 2 


while the details 


of the methodology are outlined in 


Appendix AA.3| In 



the p-spin case, one can use an alternative way to com- 
pute XFp(i) based on power counting in iV~^, the inverse 
number of spins. This provides a complementary way 
to show that dynamical fluctuations are indeed given by 
ladder diagrams. 



Let us now present our numerical results. In Fig. I12| 
we show a comparison of XFp(i) and XT{t) for various 
temperatures approaching Tc from above. Clearly, XT{i) 
is remarkably similar to Xfp(0 in this regime, exhibit- 
ing a well defined regime at short times that grows as 
a power-law with the critical mode-coupling exponent 
a — 0.395, and a well-defined power-law at later times 
that grows with the von Schwcidlcr exponent b = \. Note 
also that the height of the peak scales as t^^^ (where r 
is the relaxation time) for both functions, as predicted. 
When the transition temperature is approached from the 
non-ergodic phase, only the first regime of slow growth 
with the exponent a can be observed (not shown). These 
results represents a useful benchmark for the compari- 
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son with real liquids. Indeed, as presented in Fig. [TOl 
Xi{t) for Monte-Carlo dynamics in a binary Lennard- 
Jones mixture (where vibrational modes that may ob- 
scure the exponent a are absent) shows features strikingly 
similar to those of the p-spin model, complete with a rea- 
sonably defined regimes showing both a and b exponents 
close to Tc- 

A. 2. Exact dynamical Equations 

Following Franz and Parisi [ll|, we consider the dy- 
namic of a perturbed p = 3 spherical p-spin model evolv- 
ing with the Hamiltonian Htot{S) = H{S) - eC{S,So), 



where St is the spin state at time t, C{S,S') = 
N^^^-SiS- is the overlap function, and H{S) = 
J2i<j<k JijkSiSjSk is the unperturbed p-spin Hamilto- 
nian. The Franz-Parisi susceptibility is defined as the 
linear response of the two-point correlation function eval- 
uated in the presence of the perturbation, Ce{t,0) = 
(CiSuSo)),, as 

XFp{t) = g^' . (A.l) 

The equations of motion for Ce{t,t') and the associated 
response function GJt,t') are derived using a standard 
MSR-approach [il|,[43|. 



^ ^^^^ = - fi{t)C,it,t') + j\s f"(CS.s))G,{t,s)C,(s,t')+ I ds f'{C,{t,s))G,{t',s) 
+ (3f'{C,{t,0))C\{t' ,0) + eC,{t' ,0), 
^^^P^ = - ti{t)G,it, t') + I As f"{GS, s))GS, t') 



(A.2) 



with the damping coefficient 



A. 3. Numerical algorithm 



= T + eCS, 0) + /3/'(C,(t, 0))G,{t, 0) 

„t In the following, we elucidate the technical detail to 

/ ds {f"{Ce{t,s))G,{t,s)G,{t,s) + f'{G,{t,s))G,{t,s)} solve Eq. jXU. This is a natural generalization of an 

-'0 efhcient algorithm to solve equilibrium mode-coupling 



(A.3) 

and f{x) = We have numerically solved these 

equations using the method described below. In the limit 
of e — > 0, we retrieve the equation of motion for the sta- 
tionary state; 



dCjt) 
dt 



1 /■* 

= -TC{t) + - / ds C'^{t~s) 

2 Jo 



dC{s) 
ds 



(A.4) 



where C{t) = Ce=o(i,0). 

The temperature derivative XT{t) = dG{t)/dT is eval- 
uated by simple numerical differentiation of C{t) with 
finely spaced temperature points. 



equation developed by Fuchs et al. [49| to nonstation- 
ary systems. The meth od g iven here can be also applied 
for the aging dynamics [90|. 

First, we shall introduce a new quantity, Q^(t,t') by 



Q,{t,t') = l-C,{t,t')~ j dsG,{t,s), (A.5) 

where the subscript e has been omitted for simplifica- 
tion. This function monitors the degree of violation of 
the fluctuation-dissipation theorem. With this new func- 
tion, the MCT equation, Eq. (jA.2p can be rewritten as 



f dCe{t,t') 

dt 

dt 



fi'{t)C,{t,t') 



ds 



nt,s) 



dC,is,t') 



ds 



J l^'^) — — C^[s,t ) 



l + fi'{t)-n'{t)Q,{t,t') 



ds 



nt,s) 



dQ,{s,t') 
ds 



ds 



+ f"{t,s) 



Peit,t'), 



dQejt^s) 

ds 



{l~Q,{s,t')} 



r 



m,t') 



with ^'(t) = 1 + Pe(t, t) and 
P,it,t')=eC,it\0), 



where f{t,t') = f{C^{t,t')). In the above expression, the 
temperature T was absorbed to time, so that all quan- 



ds 



dQeit',s) 
ds 



+nt,s) 



dQe{t,s) 

ds 



C,{t\s) 
(A.7) 
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titles in the equations are dimensionless. Integration of 
Eq. (|A.6p can be implemented by discretizing the two 
dimensional plane of the times (t,t') with t > t' into a 
cubic lattice of the grid size S . Note that Eqs. (jA.6[|A.7l) 
consist of four types of time integrals; 





j ds A{t,s) 


dB{s,t') 
ds ' 




J ds A{t,s) 


dBit,s)^^ 

a C{s,t) 

OS 




/ ds A{t, s 
Jo 


dB{t\s) 
' ds ' 




/ ds A{t, s 
Jo 


.dB{t,s) 

) C{t,s 



(A., 



Likewise, other integrals can be approximated as follows. 
E ldA'^hBu-Bu-lKQ^,+Q^^,) 

■m -J 
1=3 + 1 

j 

1=1 



(A.13) 

With this discretization, the nonlinear intcgro- 
differential equation, Eq. (jA.6p . can be written in 
a form of a simultaneous nonlinear equation as 



These integrals are evaluated by discretizing the time as 
ti ~ iS and slicing into pieces as follows, /(^^(t — ti, t' = 
tj) = (i > j), for example, is written as 



V, = M,-F,(V,) + N, 



(A.14) 



z<;urds.(.,s)^ 



+ / ds A{t,,s) 



~Ai^mBm,j -^ijBjJ + ^ ^ / ds A[ti, s) 



dBjs^t^) 
ds 

dB{s,tj 



ds 



E 

i=j+i 



dAiU,s) 
ds B(s,tj), 

tl-1 



(A.9) 

where m = [(i — j)/2] is the integer closest to but smaller 
than (i — j)/2. Using an approximation. 



* ds ^4^i?(s) « {A{h) - A{t,)} X \ Tds i?(s), 
1 OS 5 Jt, 



which is exact up to 0{5'^) [4§], we arrive at 
r(i) 



(A.IO) 



where = {Cio,--- ,Cu,Qta,--- ,Qii) and Fi(Vi) = 
(/'(C.o), • • • , fiCu), /" (Qo), • • • , f"{Cu)) are {2i + 2)- 
dimensional vectors. The matrix Mj and the vector 
are functions of the friction coefficient the vectors 
at the earlier times (V;,r/) with I < i, and a set of 
bond integrals W = (dC^^), dCf^Q^'^'^, dQ^''\ df"^''\ 
df"^''\ df"'^''\ df"^"^). Equation ((XT4l) can be solved self- 
consistently using the following procedure. 

1. First, prepare the array of exact V.;, Fj, and W 
for < j < I < 7Vf/2 with a very small time grid 
5 such that NtS ^ 1 by short time expansion of 
Eq. (TOl) . 

2. For i = Nt/2 + 1 and for j very close to but 
smaller than i, we import the values of the pre- 
vious time, expecting the short time dynamics at 
{i — j)5 ^ 1 is not affected by the perturbed 
field or by aging. More specifically, we choose 
an integer TVshort <C Nt/2 and assign the values 



G, 



[h) 



ij ^ C'i-ij-i, dC^ j 
for i ~ N short < j <i- 



dCl^'^\j_i and so forth 



' — A- R . — A R 

J-ij — ^i.rn-^Jm.j 



l=m+l 

where 



{Bi, ~ Bi_,,)dA^^ - ^ (A, 

1=3 + 1 



A,i_i)dBl'f, 
(A.ll) 



3. For i = Nt/2 + 1 and for < j < i - A^short, we 
solve Eq. (jA.14p self-consistently by iteration. The 
iteration is done by choosing the initial array as 
Vi = Vi_i. The bond integrals arc calculated us- 
ing 











dA^^ 









ds A{s, tj), 



dA) 



(v) 



(A.12) 



12 

— {-A^j+2 + 8A^j+i + 5A,^j) . 



(A.15) 



ds A{U,s) 



are the integrals over the horizontal and vertical lattice 
bond, respectively (we refer to them as bond integrals). 



At every iterations of Eq. (|A.14p for V^, all ele- 
ments of the bond integrals dA^'^^'^^ and, thus M 
and N, arc updated using Eq. (|A.15p . 

4. Keep the procedure 2 and 3 for A^t/2 < i < Nf. 
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5. Once all solution for < i < Nt are obtained, we 
decimate the number of variables by half in order 
to save the memory space to explore further for the 
longer time. We discard half variables and renew 
all variables by the following rules; For V = (C, Q), 

V2^,2J^V„J. (A.16) 

For bond integrals, 

^ ^ ^ (A.17) 



Then, the time grid is doubled. 

26 6. (A.18) 

6. Repeat the procedures 2-5 with the doubled grid 
size. 

We have checked that, in order to obtain a stable result 
up to the order oit = 10^^ as shown in the present work, 
we need the number of grid of Nt = 1024 and iVghort = 32, 
starting the initial grid size oi 6 = 10^^". 
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